##################################################
## Replication Code for "Citizen Responses to Donor-Centeredness in the US-China Public Diplomacy Competition"
##################################################

#### install packages ####

packages <- c("reshape", "car","coin", "xtable","eRm", "GK2011",   "extrafont", "xtable",
              "rjags", "rlist", "maps", "dplyr", "readxl", "readr", "ggplot2", "reshape2", "gsheet",
              "tidyverse", "forcats", "countrycode", "foreign", "nnet","gsubfn", "RColorBrewer",
              "ggthemes", "cjoint","cregg", "xlsx", "stringr", "estimatr", "sjPlot","haven","jtools", "lubridate")

package_installed <-
  sapply(packages, function(pack)
    pack %in% rownames(installed.packages()))

if (any(!package_installed)) {
  sapply(packages[!package_installed], install.packages)
}

sapply(packages, require, character.only = TRUE)
rm(packages,package_installed)

#### functions ####

ate.estimator2 <- function(treat, control,
                           regression=TRUE, group=NULL, 
                           robust=TRUE, cluster=FALSE, x=NULL){
  n.treat = length(treat)
  n.control = length(control)
  if(regression){
    small.clean <- data.frame(y = c(treat, control),
                              treatment = c(ifelse(!is.na(treat), 1, NA),
                                            ifelse(!is.na(control), 0, NA)))            
    if(is.null(x)){
      ## small <- cbind(treat, control)
      
      fit <- lm(y ~ treatment, data=small.clean)
    }else{
      x <- as.matrix(x)
      small.clean <- cbind(small.clean, x)
      colnames(small.clean) <- c("y", "treatmen", colnames(x))      
      fit <- lm(y ~ ., data=small.clean)
    }
    if(robust){
      out <- lmtest::coeftest(fit, df = fit$df,
                              vcov = vcovHC(fit, type = "HC2"))
    }else if(cluster){
      out <- lmtest::coeftest(fit, vcov=vcovHC(fit, type="HC0", cluster=group))
      
    }else{
      out <- lmtest::coeftest(fit, df = fit$df) 
    }
  }else{
    est = mean(treat, na.rm=TRUE) - mean(control, na.rm=TRUE)
    se = sqrt(var(treat, na.rm=TRUE)/n.treat +
                var(control, na.rm=TRUE)/n.control)
    out = data.frame(Estimate = est,
                     SE = se,
                     t = est/se)
  }
  return(out)
}

#### data load ####

load("pd_all.Rdata")
all_data <- all_data[c("iso3c", "pd.recipient.treat.usa", "pd.recipient.control.usa", "pd.recipient.treat.china", "pd.recipient.control.china")]
all_data <- all_data[rowSums(is.na(all_data)) == 3, ] 

#### Figure 1: Experiment Results on the Effects of Donor-Centeredness: (a) US Culture Center ####

pvalue.treat12a <- all_data %>% group_by(iso3c) %>% 
  mutate(est = ate.estimator2(treat = pd.recipient.treat.usa,
                              control = pd.recipient.control.usa)['treatment',1],
         se = ate.estimator2(treat = pd.recipient.treat.usa,
                             control = pd.recipient.control.usa)['treatment',2],
         pvalue = ate.estimator2(treat = pd.recipient.treat.usa,
                                 control = pd.recipient.control.usa)['treatment', 4])%>%
  distinct(iso3c, est,se, pvalue) %>% dplyr::select(iso3c, est,se, pvalue)

all.12a <- all_data %>% ungroup %>% 
  mutate(est = ate.estimator2(treat = pd.recipient.treat.usa,
                              control = pd.recipient.control.usa, robust=FALSE, cluster=TRUE, group = "iso3c")['treatment',1],
         se = ate.estimator2(treat = pd.recipient.treat.usa,
                             control = pd.recipient.control.usa, robust=FALSE, cluster=TRUE, group = "iso3c")['treatment',2],
         pvalue = ate.estimator2(treat = pd.recipient.treat.usa,
                                 control = pd.recipient.control.usa, robust=FALSE, cluster=TRUE, group = "iso3c")['treatment', 4])%>%
  distinct(est,se, pvalue) %>% dplyr::select(est,se, pvalue) %>% mutate(iso3c = "All countries")  %>%
  relocate(iso3c)


g.12a <- ggplot(pvalue.treat12a, aes(x = est, y = reorder(iso3c, est))) +
  geom_point() +
  geom_pointrange(aes(x = est, xmin = est - 1.96 * se, xmax = est + 1.96 * se)) +
  geom_vline(xintercept = all.12a$est, alpha = 0.2, size = 13, col = "brown") +
  geom_vline(xintercept = all.12a$est - 1.96 * all.12a$se, linetype = "dotted", alpha = 0.4, col = "brown") +
  geom_vline(xintercept = all.12a$est + 1.96 * all.12a$se, linetype = "dotted", alpha = 0.4, col = "brown") +
  theme_classic() + 
  geom_vline(xintercept = 0, linetype = "dashed", color = "grey60") +
  annotate("text", x = all.12a$est, y = 14, label = "Sample Average", 
           hjust = 0.5, size = 6, col = "brown") +
  coord_cartesian(xlim = c(-0.3, 0.3), # This focuses the x-axis on the range of interest
                  clip = 'off') +   # This keeps the labels from disappearing
  labs(x = "Against US centeredness                                   For US centeredness", y = "") + 
  theme(plot.margin = unit(c(2, 3, 2, 2), "lines"),  # Adjust margin
        axis.text.x = element_text(size = rel(1.8)),  # Increase x-axis font size
        axis.text.y = element_text(size = rel(1.8)),
        axis.title.x = element_text(size = 14)) 

g.12a 
ggsave(file="dotplot_12a_all.pdf", width=8, height = 6, family="sans", dpi = 600)

## Figure 1: Experiment Results on the Effects of Donor-Centeredness: (b) Chinese Culture Center

pvalue.treat12b <- all_data %>% group_by(iso3c) %>% 
  mutate(est = ate.estimator2(treat = pd.recipient.treat.china,
                              control = pd.recipient.control.china)['treatment',1],
         se = ate.estimator2(treat = pd.recipient.treat.china,
                             control = pd.recipient.control.china)['treatment',2],
         pvalue = ate.estimator2(treat = pd.recipient.treat.china,
                                 control = pd.recipient.control.china)['treatment', 4])%>%
  distinct(iso3c, est, se, pvalue) %>% dplyr::select(iso3c, est, se, pvalue)

all.12b <- all_data %>% ungroup %>% 
  mutate(est = ate.estimator2(treat = pd.recipient.treat.china,
                              control = pd.recipient.control.china, robust=FALSE, cluster=TRUE, group = "iso3c")['treatment',1],
         se = ate.estimator2(treat = pd.recipient.treat.china,
                             control = pd.recipient.control.china, robust=FALSE, cluster=TRUE, group = "iso3c")['treatment',2],
         pvalue = ate.estimator2(treat = pd.recipient.treat.china,
                                 control = pd.recipient.control.china, robust=FALSE, cluster=TRUE, group = "iso3c")['treatment', 4])%>%
  distinct(est,se, pvalue) %>% dplyr::select(est,se, pvalue) %>% mutate(iso3c = "All countries")  %>%
  relocate(iso3c)


g.12b <- ggplot(pvalue.treat12b, aes(x = est, y = reorder(iso3c, est))) +
  geom_point() +
  geom_pointrange(aes(x=est, xmin=est-1.96*se, xmax=est+1.96*se)) +
  geom_vline(xintercept = all.12b$est, alpha = 0.2, size=12,  col="brown") +
  geom_vline(xintercept = all.12b$est-1.96*all.12b$se, linetype="dotted", alpha = 0.4, col="brown") +
  geom_vline(xintercept = all.12b$est+1.96*all.12b$se, linetype="dotted", alpha = 0.4, col="brown") +
  theme_classic() + 
  geom_vline(xintercept=0, linetype="dashed", color="grey60") +
  annotate("text", x = all.12b$est, y = 14, label = "Sample Average", 
           hjust = 0.5, size = 6, col = "brown") +
  coord_cartesian(xlim = c(-0.3, 0.3), # This focuses the x-axis on the range of interest
                  clip = 'off') +   # This keeps the labels from disappearing
  labs(x = "Against US centeredness                                   For US centeredness", y = "") + 
  theme(plot.margin = unit(c(2,3,2,2), "lines"),  # Adjust margin
        axis.text.x = element_text(size = rel(1.8)),  # Increase x-axis font size
        axis.text.y = element_text(size = rel(1.8)),
        axis.title.x = element_text(size = 14))
g.12b 
ggsave(file="dotplot_12b_all.pdf", width=8, height = 6, family="sans", dpi = 600)
